function rho = Rho(sigma, r, K, S0, dx, dt, T, scheme, type, XminFactor, ...
	XmaxFactor, dr)

[plow, vtl, Phi, Delta, Gamma, S] = PDEFD(sigma, r - dr, K, S0, dx, dt, T, ...
   scheme, type, XminFactor, XmaxFactor);

[phigh, vtu, Phi, Delta, Gamma, S] = PDEFD(sigma, r + dr, K, S0, dx, dt, T, ...
   scheme, type, XminFactor, XmaxFactor);

rho = zeros(length(vtu(:, 1)), 1);

rho(:) = (vtu(:, 1) - vtl(:, 1)) / (2 * dr);

